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Abstract 

A  numerical  study  is  made  of  the  equilibrium  properties 
and  thermal  conductivity  of  one -dimensional  lattices  with  various 
degrees  of  anharmonicity  and  coupling  defects.   The  equilibrium 
energy  flux  autocorrelation  function  is  determined,  and  its  relation- 
ship to  the  coefficient  of  thermal  conductivity  is  discussed.   The 
coefficients  of  thermal  conductivity,  together  with  other  non- 
equilibriiim  properties,  are  determined  for  several  applied  temperature 
gradients  and  lattice  conditions. 
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1.   Introduction 

Despite  the  numerous  studies  of  lattice  thermal  conductivity 
since  the  classic  investigation  of  Peierls  [l],  there  remains  a  n"amber 
of  unresolved  conceptual  and  analytic  problems  [2,3] •   The  theoretical 
clarification  of  these  problems  is  impeded  in  the  case  of  real  solids 
by  effects  resulting  from  complicated  interatomic  forces,  lattice 
structures,  impurities,  defects,  and  quantum  effects--all  of  which 
are  largely  secondary  to  many  of  the  conceptual  difficulties,  and  even 
to  many  of  the  analytic  problems.   The  purpose  of  the  present  study  is 
to  report  on  mamerical  calculations  of  the  thermal  conductivity  and 
equilibri\im  fluctuations  for  some  simple  one -dimensional  lattices  which 
have  a  minimum  number  of  complicating  factors.   Thus  our  present 
objective  is  not  to  reproduce  results  which  can  be  better  obtained  from 
real  experiments  on  real  solids,  but  to  obtain  information  about  systems 
which  are  simple  enough  for  future  detailed  theoretical  analysis.   It 
is  hoped  that  by  simplifying  the  system,  and  obtaining  detailed 
information  which  is  impossible  to  obtain  by  physical  experiments,  the 
present  results  will  aid  in  clarifying  some  of  the  conceptual' aspects 
of  lattice  thermal  conductivity  and  act  as  a  detailed  check  of  the 
standard  analytic  approximations  which  are  used  in  the  case  of  real 
solids.   It  should  be  made  clear,  however,  that  the  degree  to  which 
studies  of  the  present  type  can  aid  in  clarifying  theoretical  questions 
is  rather  limited.   The  theoretical  idealization  of  infinite  lattices 
can  only  be  roughly  approximated  by  numerical  calculations,  which  are 
at  present  limited  to  rather  small  lattices  ( 100-1000  particles)  due 
to  computation  time.   While  this  size  limitation  does  not  appear  to 


have  any  significant  effect  on  the  intensive  nature  of  the  coefficient 
of  thermal  conductivity,  it  does  have  a  profound  effect  on  the  "behavior 
of  the  energy  fl-ox  autocorrelation  function.   Because  of  this,  it  is 
very  difficult,  if  not  impossible,  to  check  the  very  interesting  and 
important  theoretical  relationships  arising  from  fluctuation- dissipation 
theories.   On  the  other  hand,  the  very  occurrence  of  these  difficulties 
indicates  certain  limitations  resulting  from  idealizations  in  the 
standard  fluctuation-dissipation  theories. 

The  model  lattice  and  reservoir  used  in  the  present  calcula- 
tions are  described  in  the  following  section.   Some  aspects  of  the 
computations  are  discussed  in  section  3-   In  section  k   a  n'omber  of 
equilibrium  properties  of  this  model  are  reported,  including  the 
behavior  of  the  energy  fl\xx:  autocorrelation  fimction.   Section  5 
contains  the  results  of  nonequilibrium  situations  involving  various 
applied  temperature  gradients  and  anharmonicities. 

2.   The  Lattice  and  Reservoir  Model 

2 

The   system  which  is   studied  is   a   linear  monatomic   lattice 

with   the  Hamiltonian 


N  N-1 

H 

i=l  i=l 


=   Y    ^  mq^  +    Y    \l  n(q.    ,-q.-£    f   -   7  V  (q.    ,-q.-^    )' 
Z_i    2        1        /__,      2   ^^1+1   ^10  3      11+1      1      o 


+  r  x- (q.   n -q.-i   ) 
h     1^^1+1  ^1     o 


(1) 


where  q.  is  the  location  of  the  i'th  particle,  and  H      is  the  equilibrium 


separation  distance.   Introducing  the  dimensionless  time  t  \l\i/in.     , 


displacements  from  equilibrium  x.  =  {^./l-    )    -  (i-l),  and  anharmonic 

coefficients  A..  =  X.    i    /\i,   k.    =  k.I   /u,  the  dimensionless  Hamiltonian, 
1    1   o'  '   1    1  o'  '  ' 

H  =  H/|i£^,  is 

N         N~l 

H  =  y  ^x^+  y  [i  (x.  -x.f   -  \   X.(x.   -x.)5  4-  I  K.(x.      ~x.)^ 
/_,    2      1    Z_.[2'i+1  i'     3   I'l+l   1     4   I'l+l   1- 

i=l       i=l 

(2) 

Typically  the  values  N  =  100,  A..  ~  10,  and  k.    ~  67  were  used  in  the 
calculations.   The  relationship  "between  these  nondimensional  units  and 
typical  physical  lattices  is  discussed  in  Appendix  A.   While  the  present 
lattice  is  monatomic,  it  can  exhibit  many  of  effects  of  "defects" 
through  the  variable  anharmonic  coefficients  X. ,    and  K. .   This  way  of 
introducing  "defects"  has  the  analytic  advantage  of  not  complicating 
the  frequency  spectrum  of  the  normal  modes,  but  it  increases  the  variety 
of  normal  mode  interactions  [h].      It,  of  course,  does  not  represent  the 
usual  defects  in  real  solids.   In  fact  the  present  method  is  intended 
to  circumvent  this  complication  while,  at  the  same  time,  probing  the 
role  of  the  anharmonic  terms  in  irreversible  processes. 

The  lattice  is  located  between  two  thermal  reservoirs  which 
are  separated  by  a  distance  (N-1)j?   or  less.   If  the  distance  is  (N-1)£  , 
as  shown  in  figure  1,  then  the  lattice  just  fits  between  the  reservoirs 
with  no  applied  pressure  at  temperature  T  =  0,   In  this  case  x   =  0  and 
X-j.  =  0  corresponds  to  contact  between  the  left  and  right  particles  and 
their  respective  reservoirs.   To  increase  the  contact  between  the 
reservoirs  and  the  lattice,  the  case  in  which  the  reservoirs  were  moved 
inward  slightly  so  that  x-,  >  0.  5  and  x  <  -0.  5  was  also  considered.   It 
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was  assumed  that  when  the  end  particles  come  into  contact  with  the  fixed 
reservoirs  they  experience  an  impulsive  collision  and  come  off  with  a 
velocity  distribution 

f(v)  =  (|v|/T)e-^  /^         (v  =  x)  (3) 

where  v  >  0  (T  -  T^ )  for  particle  1,  and  v  <  0  (T  =  T^)  for  particle  N. 
This  corresponds  to  Maxwellian  reservoirs  with  particles  of  the  same 
mass  as  the  lattice  particles.   The  factor  v  in  (3)  takes  into  account 
the  fact  that  a  slow  reservoir  particle  is  less  likely  to  "be  at  the  wall 

when  the  lattice  particle  arrives  there  (the  same  effect  as  in  molecular 

_2 
effusion).   The  mean  value  of  T  was  taken  to  he  10  ,   which  corresponds 

to  a  physical  temperature  of  roughly  300  K  (see  Appendix  A).   In  order 

to  obtain  temperature  gradients  which  are  observable  over  the  background 

statistical  fluctuations,  it  was  necessary  to  take  2(T„  -  T^  )/(T_  +  T^.  )  >  0.^. 

Therefore  the  temperature  differences  were  of  the  same  order  of  magnitude 

as  the  mean  temperature.   However,  within  the  statistical  accuracy  of 

the  present  study  (~  10'^),  no  consistent  dependency  of  the  coefficient 

of  heat  conductivity  on  the  gradient  was  observed. 

The  total  dimensionless  energy  flux  (instantaneous)  in  the 

lattice  is  given  by  [5] 

N-1  r 

J(t)  =  -  >   -  (x.  ,+x.)mx.   -X.)  -  A..(x.   -X.)   +  K.(x.   -X.)' 
^      Z_i  2  ^  1+1  ±' y    1+1  \'  1^  1+1  ±'  1^  1+1  1^ 

i-1  *- 

(^) 

This  expression  was  used  to  evaluate  the  energy  flux  autocorrelation 
function  for  the  equilibrium  system.   The  harmonic  approximation  of  J(t), 


namely 

N-1 


J^(t)  =  -  >    ^  (x.  ,-x.)(x.   -X.)  (5) 

H  '  Z_i   2  ^  1+1   1''^  1+1  i'  ^    ' 

1=1 


which  is  frequently  used  in  theoretical  investigations  [1,3^5^6],  was 

also  examined.   In  the  cases  studied,  J-rjCt)  was  not  found  to  be  the 

n         

dominant  part  of  j(t)--which  casts  considerable  doubt  on  many  theo- 
retical treatments  of  the  energy  flux.   The  average  heat  flux,  J, 
through  the  nonequilibrium  lattice  was  determined  directly  from  the 
time-averaged  energy  exchange  at  the  two  reservoirs.   This  was  used  to 
obtain  the  dimensionless  coefficient  of  heat  conductivity,  K   ,    through 
Fourier's  law 


J  =  -K^   ^T  (6) 

The  temperature  of  each  particle  in  the  lattice  is  defined  to  be  twice 

3 

the  time-averaged  kinetic  energy 


T.  =  v^(t)  (7) 


The  internal  temperature  gradient  (dimensionless)  in  (6)  is  determined 
by  a  fit  of  the  values  of  T.  to  a  straight  line.   The  internal  gradient 
is  always  considerably  less  than  the  applied  temperature  gradient, 
(T   -  T^ )/lOO,  because  of  large  temperature  differences  between  the 
lattice  ends  and  the  reservoirs.   It  should  be  also  emphasized  that  (t) 
is  simply  a  definition  of  T.,  and  does  not  imply  local  equilibrium.   It 
would  be  much  better,  but  rather  difficult,  to  relate  the  local 


temperature  to  energy  (or  velocity)  distributions.   What  this  relation- 
ship should  be  is,  however,  not  entirely  clear--the  concept  of  local 
equilibrium  may  be  too  naive  for  a  nonequilibrium  lattice. 

5.   Computational  Features 

The  problem  was  programmed  in  FORTRAN  II  except  for  the  inner 
loop  integration  of  the  dynamical  equations.   This  Runge-Kutta  integra- 
tion procedure  was  written  in  NICAP,  the  assembly  language  for  the 
University  of  Illinois  computer,  ILLIAC  II.   It  was  found  that  the  actual 
running  time  on  ILLIAC  II  was  almost  exactly  equal  to  the  dimensionless 
time  of  the  problem  itself. 

The  boundary  interactions  play  an  important  role  in  the 
selection  of  the  method  of  integration  of  the  2N  equations  of  motion 
which  derive  from  the  dimensionless  Hamiltonian  (2); 

V,  =  F.  -,  -  F.       ,       X.  =  V. 
1    1+1    1      '        11 

2  5 

where  F.  =Ax.-\.(Ax,)  +/<.('^x.)   with/!^>c.=x.  -x,  and  the  con- 
1      111      11  1    1    1-1 

dition  of  free  ends  requires  Ax,  =  Ax.,  -,  =  0.   Each  impulsive  inter- 

1      N+1  ^ 

action  with  the  reservoir  erases  the  past  history  of  the  particle 
motion  and  necessitates  a  restart  of  the  integration-   In  order  to 
allow  high  interaction  rates  with  the  attendant  better  statistical 
behavior,  the  system  is  studied  under  a  high  pressure  configuration  with 
reservoir  interactions  every  10  to  20  iterations  of  the  difference 
equations.   Rather  than  mix  integration  schemes  a  Runge-Kutta  method 
was  adopted  for  the  entire  problem  with  no  multistep  intermediate 
calculation. 


In  deciding  on  the  step  size  and  the  order  of  the  method  to 
be  employed  there  is  the  usual  compromise  between  accuracy  and  running 
time.   A  fifth  order  Runge-Kutta  method  due  to  Butcher  has  been  shown  [7] 
to  have  the  best  accuracy  of  a  typical  set  of  integration  methods  for  a 
fixed  computing  time  on  this  problem  and  in  a  suitable  error  range. 
A  step  size  of  0.1  time  uniiss  was  chosen.   It  is  by  no  means  clear  that 
this  selection  is  optimal.   The  strategy  adopted  is  a  compromise  of 
conflicting  requirements. 

On  the  one  hand,  an  important  characteristic  of  the  system 
under  study  is  the  possibility  of  performing  an  energy  check  by  balancing 
the  reservoir  energy  fluxes  against  the  calculated  energy  resident  in 
the  lattice.   The  change  in  internal  energy  will  be  exactly  accounted 
for  by  the  external  flux  except  for  errors  in  the  integration  which 
are  mainly  truncation  errors  in  our  computations.   Roundoff  error  is 
not  significant  because  of  the  IJ-significant-figure  precision  of  the 
ILLIAC  IE  computer  used  in  this  work.   Thus,  by  maintaining  high  accuracy, 
the  energy  check  can  be  used  to  detect  possible  machine  errors. 

On  the  other  hand,  if  we  relax  the  requirement  of  high 
accuracy  in  the  energy  check  and  use  a  longer  time  step  with  perhaps 
a  less  accurate  integration  scheme  the  whole  effect  of  truncation  error 
could  be  compared  with  a  physical  experiment  with  leakages  between  the 
system  and  an  external  heat  bath.   As  long  as  these  leakages  are  small 
compared  to  the  main  fluxes  in  the  system  the  steady  state  configuration 
will  not  be  sensibly  different  from  the  case  of  an  isolated  system. 
The  statistical  errors  which  arise  because  of  the  small  number  of  inter- 
actions with  the  reservoir  per  hour  of  computing  time,  about  2000  in  our 
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ccmputations,  can  be  reduced  by  increasing  the  time  step.   Since  the 
statistics  will  improve  roughly  as  the  square  root  of  the  number  of 
interactions  or  dynamical  time  a  twofold  improvement  in  statistics  will 
require  a  fourfold  increase  in  the  time  step.   Using  a  fifth  order 
method  this  means  an  error  term  three  orders  of  magnitude  larger  for  a 
halving  of  the  statistical  error.   The  advantages  of  a  small  error 
term  with  a  good  energy  check  leads  to  the  step  size  of  0.1  units  we 
selected. 

The  number  of  interactions  per  unit  of  computing  time  can  be 
increased  by  reducing  the  n\amber  of  particles.   No  improvement  is 
realized  by  this,  however,  since  the  lifetime  of  these  pulses  is  shorter 
and  we  have  fewer  particle  parameters  to  determine  temperature  gradients 
and  other  functionals.   It  would  appear  that  the  information  to  be 
obtained  is  a  function  of  the  computing  time  alone. 

A  rough  feeling  for  the  computation  involved  can  be  had  by 
observing  that  for  100  particles  we  have  200  equations  with  5  Runge- 
Kutta  steps  per  equation.   A  single  time  step  of  0.1  units  takes  about 
0.1  sec  of  computing  time  so  each  basic  integration  step  takes  100  micro- 
seconds per  particle  including  all  overhead  in  the  problem  and  in 
addition  to  the  necessary  function  evaluation  and  steps  in  the  integra- 
tion itself.   Again,  since  we  are  interested  in  statistical  quantities, 
dramatic  improvement  on  problems  of  this  type  must  await  the  advent  of 
new  parallel  computers  of  the  ILLIAC  IV  class. 


h,      Equililprium  Results 

To  facilitate  the  discussion  of  the  results  in  the  remaining 
sections,  it  is  useful  to  introduce  some  terminology  and  basic  numerical 
values: 

N  =  100;  K.    =%}?.    (Appendix  A);   T   =  ^  (T^  +  tJ  =  0.01 
'   1    3   1  av   2   R    L 

"no  defects"  means  A..  =  \  =  10 

1 

"defects"  means  that  \%   of  the  A.,  have  a  value  of  1.5^,  15^  have 
a  value  of  0.5^.,  and  ^QPJo   have  the  value  of  10  (Appendix  B) 

"pressure  applied"  refers  to  T  =  0  and  means  that  x,  >  0.5, 
=<„<-0.5 

"no  pressure  applied"  means  that  x  >  0  and  x  <  0 

"force"  (on  each  reservoir)  equals  the  time-averaged  momentum 

exchange  between  the  reservoir  and  the  end  lattice  particle 

CFL  and  CFR  is  the  collision  frequency  of  the  left  and  right 
particle  with  its  respective  reservoir 

"Applied  VT"  =  (T^  -  T  )/l00;  Internal  VT  is  obtained  from  a 
K     li 

fit  of  the  T.  to  a  straight  line. 

These  values  will  be  implied  in  all  the  following  cases  unless  other- 
wise stated.   The  time  averages  were  over  various  times,  up  to  ^0,000 
seconds  (dimensionless) ,  which  represents  ^00,000  integration  steps. 
It  should  be  noted  that  ^0,000  seconds  is  the  time  it  takes  a  sound 
wave  to  cross  the  harmonic  lattice  UOO  times  and  roughly  TOO  times  in 
the  anharmonic  lattice. 

Some  of  the  equilibrium  properties  of  the  harmonic  and 
anharmonic  lattice  were  determined,  both  for  comparison  with  the  non- 
equilibrium  cases  and  because  of  their  intrinsic  importance. 


10 


Typical  examples  of  the  equilibrium  dynamics  of  the  lattice 

_2 
between  the  thermal  reservoirs  (T_  =  T^  -  10  )  are  shown  in  Figures  2 

and  3.   These  figures  show  the  displacement  (multiplied  by  three)  of 

each  of  the  one  hundred  particles  as  a  function  of  time.   Figure  2 

illustrates  the  behavior  of  the  harmonic  lattice.   Two  features  are 

notable  in  this  figure.   First,  the  lattice  spends  a  considerable  amount 

of  time  not  in  contact  with  the  reservoirs  (i.e.,  contracted  in  length), 

and  tends  to  "slosh"  back  and  forth  between  the  two  reservoirs.   Secondly, 

the  very  energetic  disturbances  (which  trace  a  line  across  the  figure) 

have  a  velocity  very  near  the  expected  value  of  1.0  (one  lattice 

spacing/sec).   However,  when  two  of  these  disturbances  interact,  there 

is  a  delay  in  the  propagation  of  approximately  10  seconds  (indicated 

by  a  shift  in  the  heavy  lines  where  they  intersect).   Presumably  this 

delay  is  due  to  the  appearance  of  high  Fourier  components  during  the 

interaction,  with  a  corresponding  reduction  in  the  group  velocity.   The 

effect  of  this  delay  is  that  the  time  required  for  a  disturbance  to 

cross  a  harmonic  lattice  is  not  equal  to  the  length  of  the  lattice 

divided  by  the  sound  velocity,  but  somewhat  longer  (depending  on  the 

niamber  of  such  interactions). 

The  situation  in  the  case  of  an  anharmonic  lattice  with 

defects,  shown  in  Figure  J),    is  quite  different.   The  contraction  of  the 

lattice  is  sometimes  quite  large,  but  it  lasts  for  a  much  shorter 

period  of  time.   The  large  amplitude  disturbances  now  cross  the  lattice 

with  a  velocity  of  1.68  to  1.55.   Not  only  is  the  velocity  larger  than 

in  a  harmonic  lattice,  but  there  is  no  longer  any  delay  when  two  large 

amplitude  pulses  intersect.   Both  of  these  results  indicate  that  there 
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has  teen  a  considerable  change  m  the  frequency  vs,  wave  length  relation- 
ship (to  the  extent  that  this  relationship  is  significant  in  the  nonlinear 
case-'-e.go^  it  will  he  amplitude  dependent.).   The  implication  is  that 
ojCk)  is  mere  nearly  linear  m  k^  and  has  a  larger  initial  slope  than  in 
the  harmonic  case.   These  conclusions  also  agree  with  previous  theo- 
retical studies  of  these  lattices  [^j.   The  situation  is  not  simple, 
however,  for  in  the  case  of  the  anharmonic  lattice  with  no  defects  the 
velocity  of  disturhances  is  1.20  to  1.2^,  with  large  interaction  delays 
(~  lO  sec,,},  and  considerable  contraction.   The  conclusion  is  that  the 
defects  have  a  large  effect  on  smoothing  out  the  dynamics,  and  increasing 
the  velocity  of  disturbances. 

The  magnitude  of  the  statistical  fluctuations  which  occur, 
when  shorter  computations  are  used,  is  illustrated  by  the  plot  of  T. 
for  the  various  particles  (averaged  over  8,000  seconds),  shown  in 
Figure  ^,   It  can  be  seen  that  the  average  temperature  of  the  particles 
deviates  from  the  applied  temperature  by  only  a  few  per  cent.   Other 
significant  results,  which  can  be  compared  to  the  nonequilibrium  values 
in  section  5>  SiVe   tabulated  in  Table  I.   The  force  on  the  reservoirs 
fluctuates  very  ,L.ittle  for  different  runs,  whereas  the  collision 
frequencies  and  energy  flux  (which  should  be  zero  in  the  present  case) 
exhibit  modest  fluctuations  over  snorter  run  times.   These  results  are 
consistent  with  the  assumption  that  the  end  lattice  particle  has  a 
Maxwej.lian  distribution.   If  the  particle  which  strikes  the  reservoir 
has  the  distribution 
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TABLE  I 


Pressure  _  ,  Run 

Defects   applied  Force  x  10  J  x  10     CFL  CFR  time 

yes       yes         6. 71  -0.12  0.270  0.270  8,000 

yes       no          5-15  -0.75  0.210  0.201+  i+,000 
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as  it  approaches  the  reservoir,  and  the  form  (5)  as  it  leaves  the 

reservoir,  then  the  average  momentum  exchange  per  collision  is 

N/n/2  (n/t'  +  ^/t  ).   This  should  equal  (force/CF),   The  values  in  Table  I 

yield  (force/CF)  =  0.2^9     for  the  longer  run,  and  0.2^+^  to  0.251  for 

-2 
the  short  runs.   For  T'  =  T  =  10  ,  the  theoretical  value  is 

v2n  X  10   =  0.25.   Obviously  this  excellent  agreement  is  only  a  weak 

test  of  the  Maxwellian  character  of  the  distribution  function--nonethe- 

less,  it  is  of  some  interest. 

Because  of  the  theoretical  relationship  between  the  energy 

flux  autocorrelation  function  and  the  coefficient  of  heat  conductivity 

[8,9]^  considerable  effort  was  made  to  determine  this  function  to  a 

reasonable  degree  of  accuracy.   The  normalized  autocorrelation  function 

is  defined  by 


t  /  t 

o  /   ° 

A(t)  =  <J(t)j(o)>/<J^>  =  lim   /   j(t)j(t+T)dt  //    J^(t)dt 
'        t  -  00  ^o  /^o 

(8) 

where  j(t)  is  given  by  (5)-   The  numerical  approximation  used  for  (8) 
was 


11,080  / 11, 080 

A(t)  -   ^    J(5n)j(5n+T)/  ^    J^(5n)        (9) 
n:=0  /  n-0 

so  the  current  was  sampled  every  five  seconds,  for  55^^00  seconds 
(t  <  500). 

The  result  obtained  for  a  lattice  with  defects  and  applied 
pressure  is  shown  in  Figure  5»   The  first  minimum  in  Figure  5?  at 
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roughly  65  seconds,  is  very  different  from  the  corresponding  value  of 
100  seconds  in  the  case  of  a  harmonic  lattice.   In  the  latter  case  this 
first  minimum  is  related  to  the  time  required  for  a  sound  wave  to  cross 
the  system;  and  for  shorter  times  A(t)  has  the  very  simple  approximate 
form 

A(t)  ~  1  -  1.7  ^      (0  <  T  <  N)  (10) 

In  the  case  of  Figure  5^  the  first  minimum  again  corresponds  roughly  to 
the  crossing  time  in  the  anharmonic  lattice  (as  determined  from  Figure  3). 
Moreover  it  is  clear  from  Figure  5  that  the  initial  behavior  of  A(t)  is 
not  simply  a  linear  decrease  with  time,  as  in  (lO),  but  that  it  also 
contains  a  part  with  an  exponential  time  behavior.   The  hope  would  be 
that  one  could  separate  the  boundary  effect  producing  the  linear 
decrease  in  A(t)  from  the  intrinsic  relaxation-time  effect.   Geometrically 
this  would  mean  that  A(t)  should  initially  behave  in  an  exponential 
fashion  in  a  coordinate  system  which  is  suitably  rotated  with  respect 
to  the  one  in  Figiare  5-   This  simple  expectation  proves  to  be  wrong. 
Instead,  it  is  found  that  A(t)  behaves  exponentially  only  if  the 
coordinates  are  both  rotated  and  shifted  by  a  suitable  amount.   The 
significance  of  this  shift  in  the  coordinates  is  not  clear  at  present. 
In  the  transformed  coordinates  the  autocorrelation  curves  is  given  by 


A'(t')  =  .55  e"-°^^'^'      (0  <  T'  <  60)  (11) 


if  the  transformed  axis  have  the  same  scale  as  in  Figure  5«   If  the 
axis  had  only  been  rotated  (by  nearly  t^/^)  ,   then  the  coefficient  in  (ll) 
would  be  0.68  rather  than  0.55«   When  (ll)  is  expressed  in  terms  of  the 
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original  coordinates,  it  is  found  that  A(t)  satisfies  the  implicit 
equation 

a(t)  ^  0.176  -  o.oiott  +  o.82i^  exp  {-.00933t  +  1.05  (a(t)-i)3 

(12) 

(0  <  T  <  60) 

These  results  do  not  seem  to  give  any  conclusive  answer  to  the  value 
of  the  relaxation  time  for  the  energy  flux  correlation  (if  it  exists). 
The  relaxation  time  from  (ll),  namely  (.01^)   =^   70  seconds,  may  not  be 
significant  because  of  the  appreciable  shift  between  the  coordinates. 
However,  it  represents  the  only  value  which  can  be  assigned  at  present. 
The  energy  flux  autocorrelation  function  was  also  determined 
using  the  harmonic  approximation  for  the  energy  flux,  Eq.  (5).   While 
the  normalized  autocorrelation  function  bears  some  resemblence  to 
Figure  5^  the  normalization  constants  were  very  different  in  the  two 
cases.   Specifically  it  was  found  that 


11,080  11,080 

^    J^(5n)  =  160.0   ;    ^    j|(5n)  =  30.8       (lO) 
n=0  n=0 


The  large  difference  in  these  values  shows  that  the  harmonic  current 
only  represents  roughly  one  half  of  the  total  energy  flux  in  the 
anharmonic  lattice.   It  would  appear,  therefore,  that  the  frequently 
used  [1,3,5^6]  harmonic  energy  flux  bears  little  relationship  to  the 
actual  energy  fliix  in  real  lattices  that  have  large  anharmonicities. 

According  to  a  n-umber  of  theoretical  investigations  [8], 
the  coefficient  of  heat  conductivity  and  the  energy  flux  autocorrelation 
function  should  be  related  by  the  expression 
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where  the  factor  N  is  essentially  the  diraensionless  length  of  the 
lattice.   While  the  question  of  the  size  of  the  lattice  does  not  seem 
to  have  been  given  any  particular  attention,  it  is  clear  that  (ll) 
should  be  interpreted  to  hold  in  the  "thermodynamic  limit,"  N  -*  0°,  in 
order  to  insure  K     is  an  intensive  property  of  the  system.   While  this 
is  rather  trival  from  a  theoretical  point  of  view,  it  is  important  in 
the  present  context,  and  probably  in  real  solids.   The  difficulty,  in 
the  extreme,  is  illustrated  by  considering  a  finite,  but  arbitrarily 
large,  harmonic  lattice.   In  this  case  the  integral  on  the  right  side 
of  (ll)  can  be  shown  to  vanish--the  energy  flux  reverses  its  sign  within 
the  time  required  for  a  sound  wave  to  cross  the  lattice.   If  (ll)  applied 
to  this  case,  then  one  would  have  K     =   0--which  is  quite  different  from 
the  usual  statement  [l]  that  /<„  -*  °°  in  the  harmonic  limit.   Note  that 
taking  the  thermodynamic  limit  after  the  integral  does  not  change  these 
results.   On  the  other  hand,  if  one  begins  with  an  infinite  lattice 
(e.g.,  using  periodic  boundary  conditions),  then  j(t)  is  rigorously  a 
constant  for  a  harmonic  lattice,  and  the  integral  in  (ll)  is  infinite. 
It  is  clear  from  this  that  some  care  must  be  taken  in  applying  (ll)  to 
finite  lattices  (such  as  in  the  present  case,  or  in  real  experiments). 
As  can  be  seen  from  Figure  5^  there  remains  a  strong  tendency  for  the 
integral  in  (ll)  to  vanish,  even  in  the  case  of  strong  anharmonicity. 
This  tendency  is  again  a  result  of  disturbances  (or  energy)  being 
reflected  at  the  ends  of  the  lattice.   Such  a  tendency  will  persist  as 
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long  as  the  correlation  time  is  longer  than  the  time  required  for  a 
sound  wave  to  cross  the  system.   Moreover,  the  correlation  time  should 
presumably  "be  of  the  same  order  as  the  lifetime  of  sound  waves 
fphonons).   Since  sound  waves  do,  in  fact,  easily  survive  while 
crossing  real  (experimental  size)  lattices,  the  correlation  time  can 
be  expected  to  be  large  compared  to  this  crossing  time.   Thus  the 
results  shown  in  Figure  5  should  represent,  fairly  realistically,  the 
behavior  of  the  autocorrelation  function  in  real  finite  solids.   Since 
a  very  large  amount  of  cancellation  still  occurs  in  the  integral  in 
(ll),  it  seems  very  doubtful  that  the  theoretical  relationship  (ll) 
can  be  expected  to  apply  in  many  real  situations.   It  will  be  shown 
later  that  (ll)  apparently  does  not  apply  to  the  present  system. 

5.   None qui lib rium  Results 

When  the  temperatures  of  the  two  reservoirs  are  made  different 
from  one  another,  a  systematic  energy  exchange  between  the  ends  of  the 
lattice  and  their  respective  reservoirs  is  observed.   From  these  results, 
it  is  easy  to  obtain  the  time-averaged  energy  flux,  J,  through  the 
lattice.   The  more  difficult  problem  is  to  obtain  a  measurable  internal 
temperature  gradient.   The  gradient  must  be  larger  than  the  statistical 
fluctuations  (e.g.,  shown  in  Figure  ^) ,  which  either  requires  long  run 
times,  large  applied  temperature  differences,  and/or  large  anharmonic 
coefficients  with  defects.   A  number  of  cases  involving  increasing  an- 
harmonicities  were  investigated  to  .deteDrmine  the  weakest  anharmonicity 
which  would  yield  an  observable  temperature  gradient. 
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In  the  case  of  a  harmonic  lattice,  illustrated  in  Figure  6, 

the  time-averaged  kinetic  energy  of  the  lattice  particles  exhibit  a 

rather  remarkable  behavior.   Figure  6  is  for  a  run  of  ^,000  seconds, 

with  an  applied  pressure,  and  ;,T  -  'T  )  =  6067  x  10   ,   The  applied 

temperature  gradient  is  indicated  by  the  inclined  straight  line  in  the 

figure.   The  average  flux  through  the  system  is  J  =  -0.2^  x  10   ,  with 

a  force  of  .01?^+,  and  CFR  =  .06k,    CFL  =  .079  (see  Tables  I,  II  and  III 

for  comparisons).   The  most  remarkable  feature  of  Figure  6  is  the  high 

degree  of  symmetry  in  the  time-averaged  kinetic  energies  of  the 

particles--despite  nearly  60O  interactions  with  asymmetric  reservoirs. 

Moreover,  the  average  temperature  of  the  lattice  is  significantly  lower 

than  ;r  (T  +  T  ).   Other  similar  computations  indicate  that  these 
2    R     L 

features  are  real  and  not  a  statistical  accident.   At  present  there 
appears  to  be  no  theory  to  explain  these  features  of  the  none qui librium 
finite  harmonic  lattice.   The  obvious  dominance  of  a  few  modes  in 
Figure  6  may,  however,  be  due  to  the  "sloshing"  effect  noted  in  Figure  1. 
The  symmetry  is  probably  due  to  the  fact  that  pulses,  traveling  in 
either  direction,  are  transmitted  without  alteration.   Thus  the  time 
average  energy  of  all  particles  would  presumably  be  the  same  if  it  were 
not  for  this  "sloshing"  effect.   The  only  predictable  feat'uT-e  of 
Figure  6  is  that  there  is  no  'uniform  temperature  gradient  inside  the 
lattice.   Obviously  Fourier's  law  does  not  apply  to  such  a  system. 
When  a  weak  anharmonicity  (X.  =  l)  is  introduced  in  this 
lattice,  the  plot  of  T.  remained  very  similar  to  the  one  for  the 
harmonic  lattice.   The  principal  difference  is  that  the  average  tempera- 
ture of  the  lattice  is  closer  to  the  average  applied  temperature  than 
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TABLE  II 


Defects        Applied         Internal  Force  CFL  CFR  K 

VT  X   10^        VT  X   10^        -J  X   10  X    10 


yes 

6.67 

2,30 

6.23 

yes 

i+.OO 

1.27 

3.66 

no 

10.00 

2.08 

11.53 

no 

6.67 

i.i+8 

7.53 

no 

4.00 

0.853 

4.82 

6.71  0.313  0.2i+5  27.0 

6.73  0.29^+  0.254  28.8 

6.87  0.352  0.2i^i+  55.5 

6.90  0.316  0.251  51.0 

6.92  0.299  0.263  56.0 


no  6.67  2.58  8.19  8.80        OA15        0.321        51.7 
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TABLE  III 


Applied     Internal  Force     CFL      CFR      K 

5  ^  -  k  2  ^ 

VT  X  10^    VT  X  10^    -J  X  10     X  10 


10.00 

2.i+5 

7.83 

5.01 

0.259 

O.lTT 

52,0 

e.Gi 

l.i+0 

i+.50 

5.02 

0.23^ 

0.183 

32.2 
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in  Figure  6.   Also  J  =  -.55  x  10~  ,  force  =  .0262,  CFL  =  .119,  and 

CFR  =  .098.   Despite  the  anharmonicity,  no  internal  temperature  gradient 

is  observed  -  presumably  because  "the  phonon  mean-free-path  is  larger 

than  the  lattice. " 

Moderate  values  of  the  coupling  constants  (e.g.,  X.    -   5) 

produced  internal  temperature  gradients  which  are  only  marginally  above 

the  statistical  fluctuations  (for  the  available  rim  times),  so  the 

remaining  computations  were  done  for  strong  anharmonicities  {X.    =   lO). 

Actually  these  "strong'"  anharmonicities  are  reasonably  realistic  in 

terms  of  real  solids.   The  resulting  temperature  gradients  when 

T   -  T   =  6.67  X  10   are  illustrated  in  Figures  7  and  8  (both  obtained 
K     L 

frcm  24,000  second  runs).   The  steeper  straight  line  in  these  figures 
again  refers  to  the  applied  gradient.   The  only  difference  between  the 
two  cases  is  that  the  lattice  in  Figure  7  contains  defects,  whereas 
there  are  no  defects  in  the  case  of  Figure  8,   The  resulting  differences 
in  the  internal  temperature  gradient  can  be  clearly  seen.   The  increased 
coupling  between  the  normal  modes  (phonons ; ,  produced  by  the  defects, 
produces  the  expected  irreversible  effects.   In  an  attempt  to  apply 
smaller  temperature  gradients,  it  was  necessary  to  use  longer  run  times 
(40,000  seconds)  in  order  to  obtain  reliable  results.   An  example  of 
such  a  computation  is  shown  in  Figure  9-   In  this  case  T„  -T  -=  4  x  10 
and  the  lattice  has  no  defects.   While  the  internal  gradient  is  small, 
the  statistical  fluctuations  are  also  quite  small. 

The  results  from  the  above  runs,  and  several  other  cases,  are 

-2 

listed  in  Table  II.   These  are  all  for  N  =  100,  \  =  10,  T   =  10   , 

'  '   av       ' 

and  applied  pressure.   The  last  entry  in  Table  II  corresponds  to  a  case 
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with  a  weak  quartic  coupling,  which  will  be  discussed  shortly.   Since 
the  internal  temperature  gradients  in  Table  II  are  known  to  only 
approximately  ten  per  cent  accuracy,  the  resulting  values  of  the 
coefficient  of  thermal  conductivity  are  independent  of  the  gradient 
to  within  this  accuracy.   Also,  as  expected,  the  defects  lower  the  value 
of  K     considerably.   Table  II  also  shows  that  the  force  on  the  reservoirs 
is  essentially  independent  of  the  applied  gradients   On  the  other  hand, 
the  ratio  CFL/'CFR  appears  to  be  a  function  of  only  the  applied  gradient 
(see  also  Table  III),  and  not  of  the  defects  or  applied  pressure.   From 
the  collision  frequencies,  force,  and  energy  flux,  one  can  determine 
the  average  moment-um  and  energy  transferred  per  collision  at  both  the 
left  and  right  reservoirs.   An  explanation  for  these  values  is, 
unfortunately,  lacking  at  present. 

In  order  to  separate  the  cubic  and  quartic  effects,  one  case 


wa 


s  run  for  K   =  (i/^}A.  ,  rather  than  K   =  (2/3)\  .   The  results  of  this 


case  is  given  in  the  last  entry  in  Table  II.   The  value  of  K     is  con- 
siderably lowered  by  this  change,  primarily  due  to  a  large  change  in 
the  internal  gradient.   Moreover  the  force  on  the  reservoirs  is 
increased  from  . O697  to  .O88,  whereas  the  ratio  CFL/CFR  is  essentially 
unaltered.   The  present  value  for  the  quartic  coupling  constant  was 


selected  on  the  grounds  that  it  is  the  smallest  value  which  does  not 

1 
2 


12    .    -3 
lead  to  two  minima  in  the  interparticle  potential,  -x  r     -  (l/3)>^r 


+  {l/h)K   r^. 

Another  method  for  differentiating  between  the  cubic  and 
quartic  terms  is  to  allow  the  system  to  occupy  a  different  volume.   The 
results  in  Table  III  are  for  a  system  with  defects  but  no  applied 
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pressure  (at  T  =  O),  and  may  te  compared  with  the  first  two  entries  in 
Table  II.   The  increased  vol-ume  of  the  system  effectively  reduces  the 
anharmonicity,  resulting  in  an  increased  value  of  /<„•   Both  the  force 
and  the  collision  frequencies  are  decreases  "because  of  the  weaker 
contact  with  the  reservoirs. 

The  interesting  question  arises  as  to  whether  the  above 
values  of  K     can  be  shown  to  be  related  to  the  equilibriiom  energy  flux 
autocorrelation  function  by  Eq.  (ll).   If  this  relationship  is  valid, 
then  it  can  also  be  written  in  the  form 


00 

^    NT   ^0 


where  A(t)  is  the  normalized  autocorrelation  function,  (8).   Using  the 
value  in  (lO),  we  obtain  <J  >  =  l6o/ll,080.   Using  the  observed  value 
of  K     =   27  (first  entry  in  Table  II ),  (12)  would  imply  that 


/   A(t)  dT  ~  19.^  (13) 

0 

On  the  other  hand,  as  already  noted,  the  integral  of  A(t)  as  determined 

from  Figure  5  is  apparently  very  nearly  zero.   In  any  case,  over  the 

interval  0  <  t  <  500,  the  integral  of  A(t),  determined  from  Figure  5^ 

is  at  least  an  order  of  magnitude  smaller  than  required  by  (13).   It 

seems  very  unlikely  that  the  difference  can  be  accounted  for  by  the 

remaining  integral  of  A(t)  for  t  >  500.   From  these  results  it  appears 

that  the  fundamental  relationship  (ll)  does  not  hold  for  the  present 

finite  lattice. 
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Histograms  of  the  velocities  of  the  end  particles,  at  the 

time  of  impact  with  the  reservoir,  were  also  obtained.   Because  of  the 

large  statistical  fluctuations,  these  did  not  accurately  determine  the 

nature  of  the  distributions,  particularly  for  high  energies.   Thus, 

while  the  distributions  frequently  had  a  similarity  to  the  Maxwellian 

p 

-V  /2T ' 
form,  f(v)  =  (l/T')ve   '    ,  an  elementary  check  of  the  data  in 

Tables  II  and  III  indicates  that  this  is  not  exactly  the  case.   If  the 

end  particles  have  a  distribution  of  this  form,  approaching  the  reservoirs, 

the  average  momentum  exchange  with  the  reservoir  per  collision  is 

VJt/2  (  V T '  +  vT),  whereas  the  average  energy  exchange  per  collision  is 

T'  -  T.   Here,  T'  refers  to  the  end  particle  and  T  to  the  reservoir. 

These  quantities  should  respectively  equal  (force/CF)  and  (j/CF).   Using 

the  values  in  Tables  II  and  III  it  is  generally  found  that  different 

values  of  T'  must  be  used  in  order  to  satisfy  both  of  these  conditions. 

This  implies  that  the  distribution  function  for  the  end  particle  as  it 

approaches  the  reservoir  is  not  Maxwellian.   Even  when  a  single  value 

of  T •  is  found  to  suffice,  the  time-averaged  kinetic  energy  of  the  end 

particle  is  generally  not  equal  to  p  (T'  +  T),   The  values  of  this 

average  kinetic  energy  of  the  end  particles  are  (in  the  order  of  entries 

in  Tables  II  and  III) 


\ 

.85    ' 

.87 

.80 

.87 

.9^ 

.80 

.80 

.88 

^N 

1.19 

1.13 

1.21 

1.13 

1.10 

1.24 

1.19 

1.10 

No  systematic  method  of  correlating  these  various  quantities  has  yet 
been  found. 
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6.   Conclusi 


on 


The  present  study  has  not  only  obtained  a  numlDer  of  results 
which  may  he  compared  with  existing  theories,  but  it  has  also  indicated 
several  limitations  of  these  theories  and  the  present  approach.   The 
principal  limitation  of  the  present  method  arises  from  available  compu- 
tation time  (or,  what  is  the  same,  lattice  size).   While  the  present 
lattice  size  does  not  appear  to  significantly  affect  the  coefficient  of 
thermal  conductivity,  it  does  have  a  strong  influence  on  the  form  of 
the  energy  flux  autocorrelation  function.   The  comparison  of  this 
function  with  the  coefficient  of  thermal  conductivity  must  await  a  theory 
which  explicitly  takes  into  account  the  finite  size  of  the  lattice. 
Another  limitation  of  the  present  approach  concerns  the  rather  poor 
statistics,  particularly  as  regards  the  velocity  distribution  function 
of  the  lattice  particles.   It  would  be  of  considerable  interest,  for 
example,  to  know  this  distribution  for  various  locations  in  the  lattice. 
At  present  there  appears  to  be  no  adequate  theory  for  such  nonequilibrium 
distributions.   Despite  these  limitations,  sufficient  information  has 
been  obtained  for  comparison  with  present  theories,  and  possibly  as  a 
guide  for  as  yet  undeveloped  theories  (e.g.,  for  a  finite  harmonic 
lattice  between  two  thermal  reservoirs). 

We  would  like  to  express  our  appreciation  to  Richard  Hicks 
for  his  assistance  in  the  final  computations  in  this  study. 
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Appendix  A 

As  a  representative  system  take 

velocity  of   sound;   v     =   i     N/ia/in     =  3  x   10     cm/sec. 

so'  ' 

_o 

equililDri-uiii  separation;  i      =  3  x  10   cm 

-23 
mass;  m  =  30  AMU  ~  5  x  10    gm. 

Then  one  dimensionless  unit  corresponds  to 

2     2  -12 

energy:  \ill      =  mv  =  ^.5  x  10    ergs. 

time:  vm/Li  -   I    N     =   10    seconds. 
'      o'  s 

-k 
Force:  |j..£   =  1. 5  x  10   dynes 

Temperature:   mv  /k  =  3.25  x  10   K,  where  k  is  Boltzmann's 
constant. 
The  dimensional  form  of  Fourier's  law  is 

cm  -sec 
where  the  tilde  represents  a  dimensional  quantity. 
The  energy  flux  per  linear  chain  (for  simple  cubic  lattice)  is 

2   ~  ,2 


sec    (f      o         To 


so  the  corresponding  dimensionless  energy  flux  J  =  j/{\ill.    )  vfi/m 
satisfies 

J   -   -K^V  T 

~     2 /  ~  ~,   2  ,     3x 

where   K     =  K  i    /kv   .      Therefore  k  T=KT{li    /mv   j 
xxos  i_Los 

=  I  X  10-9  ~j  (_^ig^)  .  6.67  X  10-^°  K  f  (-^^£^). 

3  T       cm-sec^  T       cm-sec 

For  a   typical  value   of   K  T   =  28  watts/cm.,    this   yields    kJ£   =   .187. 
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A  pressure  of  1  atmosphere,  corresponds  to  an  average  force 
per  linear  chain  of  9  ^  10    dynes,  or  a  dimensionless  force  of  6  x  10 
Thus  a  dimensionless  force  of  .06  corresponds  to  10  atmospheres  press\ire. 


Appendix  B 

For  possible  future  theoretical  investigations,  we  list  the 
fifteen  nonlinear  coefficients  with  X..  =  15  and  K.    =  5  in  the  case  of 
defects.   They  are  respectively  i  =  ^,8, 9, 13, 23, 30, 57, ^^,51^57, 63, 6^1,73, 
91,9^  and  i  -  3, 11, 15,20,25,29,3^,36,^+7, 60, 75, 77, 82,8^4,96,98. 

Appendix  C 

The  numerical  value  of  the  normalized  autocorrelation  function 
A(5n),  for  increasing  integer  values  of  n,  is  1.000,  O.696,  0.520, 
0.390,  0.280,  0.190,  .015,  -.05^4,  -.127,  -.188,  -.250,  -.31^,  ~o356, 
-.338,  -.300,  -.250,  -.191,  -.1^9,  --.109,  --071,  -.033,  +.019,  .062, 
.097,  .1^3,  «l80,  .,202,  .213,  .20^,  .185  (for  t  -  150}. 
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Footnotes 

1.  Present  address:  Hydrospace  Research  Corporation,  Rickville,  Md. 

2.  A  recent  numerical  study  of  lattice  thermal  conductivity,  which 
investigates  the  effects  arising  from  a  disordered  anharmonic 
lattice,  has  been  made  by  D.  N.  Payton  III,  M.  Rich,  and  W,  M. 
Visscher  [h]. 

3o   It  might  be  mentioned  that  the  ratio  of  the  average  potential 
energy  to  the  average  kinetic  energy  is  approximately  0.75  when 
X.    =  10  and  T  =  .01,  so  the  time-averaged  potential  energy  differs 
significantly  from  T/2. 

k.      Numerical  values  of  A(t)  are  given  in  Appendix  C. 

5.   A  recent  study  of  a  finite  nonequilibrium  harmonic  lattice  has 
been  made  by  Z.  Rieder,  J.  L.  Lebowitz,  and  E.  Lieb  (Yeshiva 
University).   Their  system  differs  somewhat  from  the  present  system, 
and  their  results  are  not  similar  to  Figure  6. 
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